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The upwind factorizable schemes for the equations of fluid was introduced recently. 

They facilitate achieving the Textbook Multigrid Efficiency (TME) and are expected also 
to result in the solvers of unparalleled robustness. The approach itself is very general. 
Therefore, it may well become a general framework for the large-scale Computational 
Fluid Dynamics. In this paper we outline the triangular grid formulation of the factor- 
izable schemes. The derivation is based on the fact that the factorizable schemes can 
be expressed entirely using vector notation, without explicitly mentioning a particulai 
coordinate frame. We describe the resulting discrete scheme in detail and present some 
computational results verifying the basic properties of the scheme/solver. 


Introduction 

This work is a part, of effort going on at NASA 
Langley for several years towards constructing a new 
generation of the flow solvers (see, for instance T ho mas 
et alb, Roberts et al“). The key idea, that was sug- 
gested by Brandt, 3 is to use a special relaxation that 
recognizes the mixed character of a system of PDFs. 
Then each sub- fact, or of the system can be treated in 
different (optimal for it) way. It is well-known that 
t he Euler equations "consist of two different factors: 
advection and full-potential operators. The advection 
part, can be treated very efficiently, say, by the march- 
ing relaxation. The full-potential operator is of the 
elliptic type in the subsonic regime. Therefore, it can 
be treated very efficiently by multigrid. In the su- 
personic regime it becomes hyperbolic: wave equation 
with respect to the flow direction. For Mach number 
substantially larger than one, the entire system can 
be solved efficiently by marching. Nearly sonic speed 
regime can be dealt with by multigrid, but requires 
some special care. 

Solvers based on such a special (Distributive) relax- 
ation were constructed initially for incompressible flow 
and were based on the staggered grid discretizations. 
The optimal multigrid efficiency was demonstrated 
Brandt and Yavneh. 4 Staggered-grid discretizations 
exist also for the compressible flow equations. How- 
ever, they all are limited to the subsonic regime, since 
they have no shock-capturing capabilities. I he stan- 
dard shock-capturing schemes, on the other hand, are 
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not factorizable, i.e. they do not reflet*! the mixed 
character of the PDFs (see Sidilkover 5 ). Therefore, it 
is not possible to construct a relaxation of Dist ributive 
type that ran be used with these schemes. 

( -learlv. there is a need for discretizations t hat are 
both factorizable and have shock-rapturing capabili- 
ties. A factorizable upwind scheme was constructed in 
Sidilkover*’ for the case of Cartesian grids. A detailed 
description of its extension to the case of structured 
body-fitted grids is given in Sidilkover et al.‘ A set of 
numerical results was presented in Roberts et al. s 

The purpose of this paper is to present a construc- 
tion of a factorizable scheme on triangular unstruc- 
tured grids. 

Euler equations and their properties 

The non-conservative quasilinear formulation of the 
compressible Euler equations in three dimensions can 
be written using vector notation as follows 

ii • Y.s = 0 (la) 

^7- Vm + Vp= 0 (lh) 

pc 2 V ■ *7 + u • Vp = 0, (1 ( ‘) 

where .$ denotes the entropy and 

clu = dp — — 7t . (2) 

c~ 

u is the velocity vector, the pressure is given by 

P= (7 - UP W 

the speed of sound 
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In this section we recall some of t lie basic properties 
of the equations (I). It is sufficient for the purpose 
of t he analysis to assume constancy of t he coefficients. 
It is known that this system of equations is of the 
mixed type: it consists of the advection and the Full- 
Potential factors. This can be made obvious by intro- 
ducing the new set of variables. 

Recall that a vector field can be decomposed into 
solenoidal and irrotational parts 

u = V x r 4- Vo. (5) 

where o is the potent ial and 0 is the streamfunction. 
The pressure gradient is related to the gradient of the 
potential as follows 

dp = —pu • Yd$ ((j) 

Substit ut ing t he new variables 0 and 0 into the pres- 
sure equation we obtain the Full- Potential equation 

/>[(’- V-' - ( (7 • V)-]o = 0 (7) 

Note, that all the terms involving the streamfunction 
cancel out. 

Performing the variable substitution in the momen- 
tum equations gives 

pit • V (V x o) = 0 (8) 

Note, that all the terms involving the potential vari- 
able cancel out. 

Introducing a new variable vort icitv 

il - Ac (*)) 

and applying operator V x to (8), we obtain 

pit YD = 0 (10) 

This verifies indeed that the Ruler system is of the 
mixed type. The advection factor is represented by 
the equations for entropy (la) and vorticity (10). The 
full-potential factor is given by (7). It also makes it 
clear t hat (for the linear constant coefficients case) the 
momentum equations ( lb) drive the solenoidal part of 
the solution, while the irrotational part of the solution 
is subject solely to the pressure equation (lc). 

In a general nonlinear case (away from singularities, 
like shocks and contact discontinuities) there is a weak 
coupling between different factors due to the so-called 
subprincipal terms. T his coupling can be neglected 
for t he purpose of t he construction of a fast solver (see 
Brandt 3 ). Therefore, the latter can rely entirely on 
the analysis of t he linear case. 

Preparations for the scheme 
construction 

When constructing a discrete approximation to the 
Filler equations, the central scheme can serve a ba- 
sic building block. However, it is crucial to include 


a certain artificial dissipation in the discretization for 
stability reasons. One of the additional problems than 
becomes how to compensate for t he loss of accuracy 
due to the artificial dissipation. Various ways to deal 
wit.li this issue received an extensive coverage in the 
literature. ( 'oust ructing a factorizabh scheme implies 
resolving this issue in a very specific way. 

FDA analysis 

The First Differential Approximation (FDA) (or the 
modified equations) corresponding to a certain discrete 
scheme is the PDFs augmented by the leading error 
terms. 

We shall start our analysis with formulating the 
FDA for the factorizable genuinely multidimensional 
scheme. The observation made in regarding the gen- 
uinely multidimensional upwind scheme introduced in 
Sidilkover was that a part of the artificial dissipa- 
tion present in the discrete momentum equations (in 
subsonic case) is proportional to the gradient of the 
residual of the pressure equation. The artificial dissi- 
pation of the pressure equation is proportional to t he 
divergence of the residuals of the momentum equa- 
tions. A vector formulation of t he entire scheme (on 
the Cartesian grids) is given in Sidilkover. 10 

The fact that the entire scheme can be expressed 
using the vector notation appeared to be very instru- 
mental for the purpose of extending the factorizable 
to the structured body-fitted grids (See Sidilkover et. 
aP and Roberts et a I s ). It is of very important for the 
purpose of this paper too. 

Tl ie FDA of a factorizable scheme for the Euler 
equations is given by the following 

q.s = 0 (11a) 

pqil -P Yp — - Y(pc 2 Y ■ it + it • Vp) — YD 

2 c 

= 0 (lib) 

pc 2 V • i7 + i7 Y]> 7 ^-rV ■ ( pit Yu + V/>) 

= 0 (Mr) 

The term YD in the momentum equations plays an 
important role when the operator q is discretized us- 
ing an advection scheme of a certain type to maintain 
the second order accuracy and factorizability of the 
whole Euler scheme. This special type of the advection 
scheme allows to upgrade the accuracy of the advec- 
tion factor to the second order without affecting the 
discrete full-potential part . For now we omit the term 
V D and consider q to correspond to the standard first 
order accurate advection scheme. 

A factorizable scheme corresponding to the FDA as 
given by (11) is stable for subsonic rase only. The 
scheme can be extended so it will be valid for tran- 
sonic/supersonic regime by a simple modification: in- 
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troduct ion of a cut-off parameter on the terms involv- 
ing pressure in the artificial dissipation 

q.s - 0( 12a) 

fiqS + Y P - ^ - Y(pc 2 Y • u + - U ■ Y P ) -YD 

2 C K 

= 0( 12b) 

pc 2 V • *7 -|- it • Vj> j-cV (/'** * + K'Vp) 

= 0,(1 2c) 

where 

k = max(l , A/ 2 ) (1*3) 

We shall demonstrate the factorizability property of 
FDA corresponding to the new scheme. Introduce the 
auxiliary variables - potential o and streamfunct ion V' 
and substitute the following expressions for the vari- 
ables it , p 

(T p l 1 _ _ 

it = (I -u ■ V)V x f 

2 c 

+V(/- (Ma) 
2 c 

dp — —p(it • V ^-rV J )dcj ( 14b) 

into (11). Introducing the vorticity variable 

O.Vx(l^-fV)Vxr (15) 

2 r 

and applying t he V x operator to the momentum equa- 
tions we obtain 

pit ■ VQ = 0. (10) 

Note, that, as well as in the PDE case, all the terms 
involving the potential variable canceled out. Substi- 
tuting the auxiliary variables into the pressure equa- 
tion, it is easy to verify that all the terms involving 
the streamfunction cancel. 

We can summarize that due to the specific form of 
the artificial dissipation, the FDA (11) of the discrete 
scheme is factorizable - it reflects the mixed character 
of the original system of PDE^s. 

Special discrete operators 

We consider the two-dimensional case from now on. 
When discretizing the derivatives, a special care needs 
to be taken of what kind of discrete operators are used 
in order to preserve the factorizability property at. the 
discrete level. 

Note, that when demonstrating the factorizability 
property of the scheme's FDA (11) we used the facts 
of the following type 

(\r rOy ~~ Oj-yOjc 

Oyy O.r Ojcy Oy ( 1 * ) 

OjryOjry = OjrxOyy 


In order to obtain the factorizable discrete scheme, 
we need to introduce some finite differences that pos- 
sess the property analogous to (17). Such finite differ- 
ences were int roduced for the case of structured grids 
in.*’ 

We have to find such differences for the case of 
triangular grids (see Figl), where (j\t/) is a local 
(non-orthogonal) frame. We shall illustrate this on 
a simple example. Consider approximating the par- 
tial derivative 0 rr . Assuming that we have at our 
disposal the "compact'' stencil that involves 7 points: 
0. 1.2, 3. 4, 5,(3. there is only one wav of doing it, 
namely by defined as follows 

0f V it — (im - 2 «o + u\)/lr (18) 

Differences defined in such a way do not have the 
property (17). We can conclude that using the com- 
pact stencil 7-point stencil only there is no wav to 
achieve this property. We know from 1 ' that the 9-point 
box structured grid stencil is sufficient for this pur- 
pose. There are several wavs to augment the compact 
7-point, stencil to the 9-point one in the current trian- 
gular grid context. It is clear that the 7-point stencil, 
therefore, needs to be augmented. We can do it by 
adding to it () more nodes: 7.8,9, 10. 1 1. 12. I he par- 
tial derivative c) jr ran then be approximated by a irtdt 
difference 

0j. v = [(w4 — 2 mo 4- mi )/2 

+ (Mio + «3 - Mv — u«)/8 (19) 

T ( m. 1 1 T U(] — Mr> — u 7 )/$]/h 

The derivatives c) yy , 0 r , 0 f/ can be approximated in the 
analogous way. It is easy to verify that such differences 
possess the property (17). 

Structure of some artificial dissipation terms 

The central part of the scheme is the constructed 
in the standard way. The artificial dissipation terms 
corresponding to the advect.ion scheme are evaluated 
in the standard fashion as well. A special care needs 
to be t aken of the other artificial dissipation terms. 

Recall, that the artificial dissipation terms in the 
momentum equations (11b) that are subject to the 
gradient operator are the residual of the pressure equa- 
tion. Denote them 

H r = pc-Y ■ it + m * V p (20) 

and the expression subject to the action of the diver- 
gence operator in the pressure equation (11c) is the 
residual of the momentum equations 

fl m — pH, ■ Vi7 + V p (21) 


3 OF 7 


American Institute of Aeronautics and Astronautics Paper 2001-2575 


9 



Fig. 1 Computational grid segment and a control volume. 


Constructing the discrete scheme 

Our goal now is to derive a discretization of a con- 
servation law (scalar and a system). For this purpose 
we need to evaluate the numerical fluxes through each 
of the faces of the dual-median cell (see Fig.l). 

Global and local coordinate frames 

11 = ( 4 4 ) ' (22) 

where and (o r/ ,/^) are the unit vectors in 

the direction of the £ and // coordinate axes respec- 
tively. The relationship between the Cartesian and 
rontravariant velocity components is described as fol- 
lows 



The .Jacobian of this coordinate rotation 



6 


Fig. 2 Nodes used for evaluating the flux through 
a face of dual- median cell. 


J — (let // — (25) 
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The inverse of the Hessian H 

h- ] = ~ 

J \ ->h 


C26) 


It is convenient to nse the sealed cont ravariant velocity 
romponents 


l — .] It V-djj CO yj 

r = jv = —uk + ro f 


(27) 


Tli relationship between the ('artesian and covariant 
velocity components is given by the following 

,7 / « \ _ I & 


IV 


V 


(28) 


U = u<\$ + 

V — no ,j + vil,) 


follows 


IV II 


jjT H _ I n ( + :j l <v £ n 'i + l) Vk 


OfO,, + o-; + i' 

The total velocity squared 

| «|“ — rr -f c 2 = ti// 4- i’V 
Scalar advection 

Consider a scalar advection equation 
T US .r + = 0. 


(32) 


(33) 


The discrete equation to solve for * at point 0 is ob- 
tained by balancing fluxes through the surface of the 
dual median cell (see Fig. 1 ) 


A\i 

£[/ /*,].- = o, 

i = l 


(34) 


the numerical flux 


/ - + M r (- s o + *2, 


where stands for a divided difference 
h^djrS — S n — •*>() 


(36) 


Euler system 

Integrating the Euler equations in the conservative 
form over the control- volume (dual- median cell) and 
applying the Green's theorem, we obtain 


F • ndl = Lb 


(37) 


vc 


where OC is t he control volume's boundary n is a unit 
vector normal to the boundary and 

\ ( pr \ 

i + 


/ 


F = 


fm 

fnr + p 


pv 

pin’ 

pv- + p 

\ (£ + ;>)«■ } 


J 


(38) 


The discrete equation to solve for s at point 0 is ob- 
tained by substituting numerical fluxes evaluated by 
analogy to (35) on all the faces of the dual-median cell 
and substituting them into the flux-balance equation 
(34). 


pur 

\ (E + p)u ) 

A conservative discretization of the Euler system ('an 
be written in the following form 

A Ki 


(29) 

t— 1 

(39) 

related as 

The numerical flux through a face can be represented 
as a sum of central and the artificial dissipation parts 

(30) 

F = F c + F d 

The central portion of the flux is given by 

(40) 

) (31) 

F e = i[F(u L ) + F(u w )]m 

(41) 


the face. The rest of this section is dedicated to the 
question of deriving the diffusive portion of the numer- 
ical fluxes. 

We would like to emphasize that it is fairly simple 
to implement, the new discretization within the exist- 
ing control- volume computer codes. It requires only 
the new numerical flux routine. Such a routine can be 
written in several simple steps, starting from the stan- 
dard upwind scheme and performing the modifications 
gradually. 

Tin standard upwind scfn nn 

The first step is to rewrite the standard upwind for 
the subsonic case without explicit mention of the char- 
acteristic variables, etc. The scheme is given by the 
following numerical flux 


(42) 


The derivation presumes that these fluxes correspond 
to the local orthogonal coordinate frame associated 
with the cell-face. Therefore, the momentum equa- 
tions diffusive fluxes can be rotated to the global co- 
ordinate frame (.r,;*/) bv in the following way 



/ /. \ 


\u\d^ \ 

■ r >) 

Ms 

= -5"< 

pcV’r + r/cd'ip 
P \U\cp' 


V h ) 


\ pcUdj! V + cd£p j 


h 

h 


h 

h 


( 43 ) 
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The obtain diffusive fluxes correspond to the quasi 
linear nonconservative formulation of the Euler equa- 
tions. The need to be transformed, therefore, to the 
form appropriate for the conservative discretization 


p'< = up 1 


(44) 


where 


1 U 0 

u 1 0 

r 0 1 

( it ■ u)/ '2 u r 


i/r- \ 

(//(■- 
r/c- 

l/h - l) + (t7-ff)/(2f 2 ) / 
(45) 


A modified, scheme 

As an intermediate step towards constructing tin* 
factorizahle scheme we can consider the following case: 



h \ 


\e\op \ 

f d = 

n 

= -\ b < 

ped^ v + v/cd'iv 

p\U\<)'{V 


/> J 


per d h ( r + (■()'• p / 


The main difference from the standard scheme is that 
the momentum equations diffusive fluxes (the second 
and third components) are now attributed to the mo- 
mentum equations in the covariant directions (£ and 
//). Therefore, the transformation back to the global 
ort hogonal frame takes the following form 


h 

h 



(47) 


This change is necessary towards eventually obtaining 
the approximation for the gradient operator term in 
(lib) 

Anotfu r nit< rmt rfiate step 



/, \ 


(C| op \ 

f J = 

h 

h 


pedp' + 1 ' / cd£ p 

f>\e\o^v 


/> / 


pev op- + p-op, / 


The coefficient in front of the pressure difference in 
the fourth component constitutes the change from the 
previous case. It is necessary in order arrive later to 
the approximation of the divergence term 

Tin factonzable scheme 

Now we shall incorporate the correction (mixed 
derivative) terms into the scheme, so that some of 
the terms in the artificial dissipation can be viewed as 
residuals of the moment um and the pressure equations 
and, therefore, are second order small. The latter mod- 
ification together with int roducing the wide differences 
that possess the commutativity properties results in a 


factonzabh scheme. Implementat ion of all these steps 
can be done in several steps as well, testing the routine 
after each modification. 

We start from re-interpreting of some standard no- 
tions. A standard "narrow" divided difference can be 
defined also as 


on 1 2 ill) 1 2 , CO 2 ;j ; )0 23 

n h ~ * G(i + 

f £012 q_ £023 


(49) 


Define a "wide divided difference" 

d f { = (25 01L 7T mi ' + 

+S ()C)l dh [ -f 
+s- [ *c)f 3 + s" 34 df 34 )/ 

(2,S' ni - + -f 4 v ),:i + .s' 18 - -f S~ 9:i + .V 03 * 1 ) 

(50) 


Similarly, we can define and . Adding the some 
specific /^-derivative terms (all the differences used here 
are narroiv) to the artificial dissipation of the pressure 
equation obtain the residual of the momentum equa- 
tion in the direction normal to the cell face 

r{\ = p{vc^ + va£)0 + + (o e n, + .yf.yo^jp 


We also add some //-derivative terms to the artificial 
dissipation of the momentum equation in £ direction 
to obtain the residual of the pressure equation. All 
the differences used here are wide. The notation (IR 
for the residual) reflects this fact 


K = (K-^V + tfV) + (Id 1 ’ + vty)p (52) 

The artificial dissipation then takes the following form 


/' = 


h \ 

f2 

h 


\_ 

2 


h tie Id's* \ 

rr,„l + h(p\0\d%M] 

hptf 1 l^v 

/ 

(53) 


1 he need for rescaling the artificial dissipation in order 
to avoid the quasi-ellipticity of the full-potential fac- 
tors approximation for the low-speed flow regions was 
established in/ 5 Tin* scaling parameters < 7 m ,<r P serve 
this purpose. I 11 the conducted preliminary numerical 
test (see next section) they were taken to be i. The 
parameter / was taken to be equal to h^. 

There remains a need to introduce some modifica- 
tions into the central part of the scheme. It is neces- 
sary for factorizabihty that the pressure gradient term 
in the momentum equations is approximated by dif- 
ferences. Using wide differences to approximate the 
pressure advection operator in the pressure equation 
is not necessary for factorizabihty but is still benefi- 
cial since it results in a better form of the discrete 
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Fig. 3 The computational grid. 



Fig. 4 Solution: density contours. 


full- potential factor. These modifications can be in- 
troduced in a very simple wav using a trick by Tom 
Roberts: adding certain terms to the artificial dissipa- 
tion. Introduce the following undivided difference 


[(2 5 02: V)f 3 + 

-(2 s (]l -c)" Y - + s Q « [ dy n 

2 lH 


+* i8a oi/ 


(54) 


The artificial dissipation terms are then augmented 
as follows 


h 

u 


f'l + v V 

h + ) 


(55) 


Returning to the global (Cartesian coordinate frame 

( /: ) ==<« 7) " ( /0 

and converting to the conservative form 

F d - AIf d (57) 


This describes the scheme that was used in the pre- 
liminary numerical experiment reported below. 


Preliminary numerical results 

The work on implementing the new scheme within 
the FUN2D 11 code has just began. Our very first aim 
is just to implement the numerical fluxes and to verify 
the correctness of the residual evaluation. 

The testca.se presented is a subsonic flow (Mach = 
.2) in a channel with a bump. The grid consists of 
1375 nodes (see Fig. 3). A second order version of the 
new scheme was used. The contour plots of density 
are presented in Fig. 4. We also present for compar- 
ison in Fig. 5 density contours of the solution to the 
same problem using the standard second order scheme 
with limiters switched off. order upwind scheme. The 
solution obtained using the new scheme is at least as 
accurate as the one using the standard scheme. 


Fig. 5 Solution obtained using the second order 
upwind scheme: density contours. 
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